Verwey-Type Charge Ordering and Site-Selective Mott Transition in Fe4O5 under Pressure

The metal–insulator transition driven by electronic correlations is one of the most fundamental concepts in condensed matter. In mixed-valence compounds, this transition is often accompanied by charge ordering (CO), resulting in the emergence of complex phases and unusual behaviors. The famous example is the archetypal mixed-valence mineral magnetite, Fe3O4, exhibiting a complex charge-ordering below the Verwey transition, whose nature has been a subject of long-time debates. In our study, using high-resolution X-ray diffraction supplemented by resistance measurements and DFT+DMFT calculations, the electronic, magnetic, and structural properties of recently synthesized mixed-valence Fe4O5 are investigated under pressure to ∼100 GPa. Our calculations, consistent with experiment, reveal that at ambient conditions Fe4O5 is a narrow-gap insulator characterized by the original Verwey-type CO. Under pressure Fe4O5 undergoes a series of electronic and magnetic-state transitions with an unusual compressional behavior above ∼50 GPa. A site-dependent collapse of local magnetic moments is followed by the site-selective insulator-to-metal transition at ∼84 GPa, occurring at the octahedral Fe sites. This phase transition is accompanied by a 2+ to 3+ valence change of the prismatic Fe ions and collapse of CO. We provide a microscopic explanation of the complex charge ordering in Fe4O5 which “unifies” it with the behavior of two archetypal examples of charge- or bond-ordered materials, magnetite and rare-earth nickelates (RNiO3). We find that at low temperatures the Verwey-type CO competes with the “trimeron”/“dimeron” charge ordered states, allowing for pressure/temperature tuning of charge ordering. Summing up the available data, we present the pressure–temperature phase diagram of Fe4O5.


Experimental Methods
Samples of Fe 4 O 5 were synthesized using a 1200-tonne Multi-Anvil Press at BGI [1] at HP-HT conditions from stoichiometric mixtures of fine powders of Fe 3 O 4 and iron (both from Aldrich company, 99.999% purity). These syntheses were performed at pressures of about 14 GPa over several hours. High-quality single crystals with linear sizes of about 20-250 µm were synthesized at high temperatures in the range of 1350-1450 • C [2]. We employed a standard multi-anvil assembly that included a Re cylindrical sample capsule, a LaCrO 3 heater, W3Re/W25Re thermocouple, and other components packed inside an octahedron made of 5% Cr 2 O 3 -doped MgO [3]. The synthesis procedure was similar to the one described in previous works [2,4]. The chemical composition of the samples was verified in a microprobe analysis using a JEOL JXA-8200 electron microscope. The crystal structure of the samples was determined using single crystal and powder X-ray diffraction (XRD) on a high brilliance Rigaku diffractometer (Mo K α radiation, λ = 0.7108Å).

X-ray diffraction studies
Data collection Run 1. The data collection was carried out at the ID15B beamline at the European Synchrotron Radiation Facility (ESRF), Grenoble, France (MAR555 flat panel detector, λ = 0.41114Å, beam size 10(V) × 10(H) µm 2 ). Sample-to-detector distance, coordinates of the beam center, tilt angle and tilt plane rotation angle of the detector images were calibrated using Si powder. XRD wide images were collected during continuous rotation of DACs typically from -20 • to +20 • on ω; while XRD single-crystal data collection experiments were performed by narrow 0.5 • scanning of the ±32 • ω range. A black single crystal of Fe 4 O 5 with dimensions of 0.02 × 0.02 × 0.01 mm 3 together with a small ruby chip (for pressure estimation) were loaded into LeToullec membrane-driven type DAC [5] equipped with Boehler-Almax diamonds with 300 µm culet size. A hole with diameter about 150 µm in steel gasket pre-indented to 40 µm was served as a pressure chamber. Neon was used both as a pressure transmitting medium and as a pressure standard [6].
The DAC was compressed to 46.8 (5) GPa with wide and step scan images being collected at each pressure point. Processing of XRD data (the unit cell determination, integration of the reflection intensities, empirical absorption correction) were performed using CrysAlisPro software [7]. A single crystal of an orthoenstatite [(Mg 1.93 ,Fe 0.06 )(Si 1.93 ,Al 0.06 )O 6 , space group P bca with a = 18.2391(3), b = 8.8117 (2), and c = 5.18320(10)Å], was used to calibrate the instrument model of CrysAlisPro (sample-to-detector distance, the detector's origin, offsets of the goniometer angles and rotation of the X-ray beam and the detector around the instrument axis). This run was carried out at the beamline 13-IDD of the Advanced Photon Source (APS), Argonne National Laboratory (Argonne, IL), with a wavelength of λ = 0.2952Å, and a spotsize of 3(V)×4(H) µm 2 . Detector tilts and sample to detector distance were calibrated using a LaB 6 standard. The laser heating process formed a few single-crystal grains along with some powder of Fe 4 O 5 . Then, we compressed these samples up to ∼100 GPa and collected both single-crystal (ω rotation of −33 • to +36 • with 0.5 • steps) and powder diffraction data.
Single-crystal diffraction processing was performed in the same manner described in Run 1, and using an orthoenstatite standard from the same batch.
We employed DACs with diamond anvil culet sizes of 300 and 150 µm, respectively.
Together with the sample in the same cell we loaded a ruby ball (for the first run) and a chip of gold (for the second run) that was used for pressure determination [6]. The DACs were filled with a Ne pressure-transmitting medium. At each pressure point we acquired the single-crystal X-ray diffraction images upon a continuous rotation of the cell with the sample around the vertical ω-axis with a step of ∆ω = 0.5 • and an exposure time of t = 0.5 s/frame.
The diffraction data were collected by a Perkin Elmer XRD1621 detector at DESY and a Pilatus 1M CdTe detector at APS. We analyzed these data with CrysAlisPro software [7] and solved the crystal structures using JANA2006 software [8]. We note that the single-crystal data were obtained up to ∼94 GPa and the images show that under pressure the structure does not change (see Figs. S2 and S3). However, at elevated pressures, particularly above 50 GPa, the quality of the SC data deteriorates. This could be due to two main reasons: 1) The crystal experiences more stress as the pressure medium becomes less hydrostatic. 2) Electronic or structural transitions often induce such lowering of the crystal quality [9], with the quality possibly improving upon completion of the transition, as seen in this case at 93.9 (14) GPa [Figs. S2 and S3]. Correspondingly the error estimates of the atomic positions and unit-cell parameters obtained above 50 GPa are large; accordingly we do not plot them.
The powder diffraction data were first integrated using DIOPTAS [10] and then analyzed by LeBail extraction using GSAS II [11]. Powder XRD was used mainly for obtaining the lattice parameters and equations of state (EOS) at the whole pressure range. The pressure uncertainties are ∼0.3-1.5 GPa from low to high pressures. Due to the large number of overlapping peaks, the c axis values were fixed according to the position of the (002) peak, which is well separated from other peaks. The reported uncertainties are given according to the standard errors obtained from the respective software used for fitting the data.

Structure solution and refinement of Fe 4 O 5
The crystal structure of Fe 4 O 5 reported by Lavina et al. [12] (space group Cmcm, Z = 4) was used as a starting model for the refinement. Since a body of the diamond anvil cell shadows more than 50% of the diffraction reflections, the reflection datasets were incomplete.
In order to improve the data/parameter ratio, we refined atomic thermal parameters in isotropic approximation. The structures were refined by full-matrix least squares against F 2 using the SHELXL-2014/7 [13] implemented in OLEX2 software [14]. The resulting R 1 values varied from 6.0 to 7.8% The most intense residual electron density peaks were not exceeding 1.8 e/Å 3 and were found within 1Å of the Fe-atoms. The detailed summary of the crystal structure refinements along with information on unit cell parameters, atomic coordinates and isotropic displacement parameters are summarized in Tables S3-S5. To evaluate the oxidation states of different Fe ions in the crystal structure we applied a bond-valence-sum (BVS) approach [15] based on the analysis of Fe-O bond lengths around the Fe ions. In this method the oxidation state of an ion V i is the sum of all its bond valences S ij , each of which is determined as S ij = exp((R ij -d ij )/b 0 ), where d ij is the distance between atoms i and j, R ij is the empirically determined distance for this cation-anion pair and b 0 is an empirical parameter, which is normally about 0.37Å [15].

Electrical transport measurements
Electrical transport measurements up to 100 GPa were performed using TAU pistoncylinder DACs. A crystal with dimensions of about 70 × 70 × 20 µm was loaded into a cavity of 100 µm in diameter, drilled in a rhenium gasket. The gasket was insulated with a layer of Al 2 O 3 -NaCl (in the atomic ratio of 3:1) that was mixed with an epoxy. Some amount of a CaSO 4 powder was prepressed into the cavity before the crystal loading to ensure a quasi-hydrostatic environment around the sample. Six platinum triangles, which served as electrodes, were placed on the culet. This assembly permitted the measurements in various DC four-probe arrangements at a given pressure. The Pt electrodes were connected to exterior conducting wires with a silver epoxy. The electrical resistance was measured as a function of pressure and temperature (for both compression and decompression cycles) using the standard four-probe method in a custom-made cryostat. At each temperature point, the voltage was measured as a function of a series of applied electrical currents to determine the electrical resistance from a slope of their dependence. The temperature was measured using a Lakeshore Si (DT-421-HR) diode in proximity to the DAC. A few ruby pieces were placed in the central region of the culet between the Pt electrodes around the sample for determination of the pressure. The pressure inside the DAC was determined both before and after each measurement using the ruby R 1 fluorescence line as a pressure marker together with the calibration scales reported in Ref. [16]. The pressure gradients in the DAC were found to be rather insignificant. They amounted to ∼5% for the distances across which the voltage was measured (15-20 µm between the tips of the Pt electrodes).

X-ray diffraction results
Our XRD experiments find the conservation of the original orthorhombic CaFe 3 O 5 -type structure up to ∼100 GPa (see Figs. S2-S5) but reveal a number of features in its compression behavior. Similar to [17], we find a steeper decrease of the unit-cell volume with pressure increase above ∼50 GPa, with ∆V /V ∼ 2.4%. Above 56 GPa the volume variation becomes more gradual and shows only an additional tiny volume drop of ∼0.3% at about 84 GPa.
Analyzing the V (P ) behavior we can highlight the existence of three pressure ranges, namely, from ambient pressure to about 50 GPa (labelled as LP phase), 56-81 GPa (HP1 phase), 6 and 84-100 GPa (HP2 phase) (see Table S1). As we show in the main text these pressure ranges correspond to different electronic states within the same crystal structure. The unit-cell volume data for these ranges can be fitted using the third-or second-order Birch-Murnaghan (BM3/BM2) equations of state (EOS) [18] (see Table S2). The fit of the V (P ) data was performed with the software EosFit GUI [19], using weights for both P and V (see Table S2).
We note a significant "softening" of the material beyond ∼50 GPa. Thus, the bulk modulus of Fe 4 O 5 almost halved and V 0 increase by ∼8%, compared to the LP phase (see Table S2). Such a behavior contradicts with the anticipated hardening of the material after the volume collapse. Similar behavior was reported for PrFeO 3 [20], in which a transition in Fe ions to the LS state was not completed simultaneously with a sharp isostructural phase transition around 50 GPa and continued as a sluggish second-order HS-LS transition at higher pressures.

Pressure-induced high axial anisotropy of Fe 4 O 5
We note an unusual compression behavior of the parameter b leading to a high axial compression anisotropy of the LP phase (see Fig. S6). In the previous work [2], it was documented that for a single crystal of Fe 4 O 5 the parameter b has a negative thermal expansion.
This anomalous behavior was tentatively addressed to an antiferromagnetic spin alignment along this axis, which can gain in the exchange energy upon a Fe-Fe distances increase [17,21]. Such spin alignment may explain also the observed unusual axial anisotropy of the LP phase and indications of the reduced compressibility of the FeO6 octahedra comprising of the Fe2 sites. We note that the high axial anisotropy collapses concurrently with the partial spin transition at ∼55 GPa, tentatively as the result of a significant decrease of T N and H hf accompanying usually HS-LS transitions [22][23][24][25][26][27]. are semiconducting-like, exhibiting a negative R(T ) slope (Fig. 2a). The curve at different fixed pressures for certain temperature ranges can be well described by an expression as follows lnR = lnR 0 + E a /k B T , where k B is the Boltzmann's constant, E a the electrical transport activation energy [see Fig. 3

(a)].
We note that a high-pressure low-temperature phase diagram of Fe 4 O 5 determined in a previous study [4] established the existence of novel charge-ordered phases at moderate high pressures up to 40 GPa. Hence, the temperature dependencies of the electrical resistance in this pressure range should correspond to a mixture of the room temperature and lowtemperature partially charge-ordered phases. The lnR(1/T ) dependences for pressures 4.

Theoretical Methods
We supplement our single-crystal (SC) and powder (PWD) synchrotron x-ray diffraction (XRD) and resistance measurements by the density functional theory+dynamical mean-field theory (DFT+DMFT) electronic structure calculations [28,29]. We study the electronic structure, magnetic properties, and valence configurations of iron of paramagnetic Fe 4 O 5 using the fully charge self-consistent DFT+DMFT method implemented within plane-wave pseudopotentials [25,26]. By using DMFT it becomes possible to treat the effects of electronic correlations (e.g., local spin and charge fluctuations) on the electronic properties of strongly correlated electron materials near the Mott transition [22,23,25,26,[28][29][30].
In our DFT+DMFT calculations we construct a basis set of atomic-centered symmetryconstrained Wannier functions for the partially occupied Fe 3d bands evaluated within DFT [31]. The DFT+DMFT calculations are performed in the local basis set determined by diagonalization of the corresponding Fe 3d occupation matrices. We used the continuous-time hybridization-expansion quantum Monte-Carlo (segment) algorithm to solve the realistic many-body problem [32]. The DFT+DMFT calculations are performed in the paramagnetic state at an electronic temperature of T el ∼ 770 K (our calculations at 580 K give quanti-8 tively similar results). In DMFT we use the Hubbard interaction U = 6 eV and Hund's exchange coupling J = 0.89 eV for the Fe 3d states in accordance with previous estimates [22,23,25,26]. The U and J values are assumed to remain constant upon variation of the lattice. The Coulomb interaction was treated in the density-density approximation, spin-orbit coupling was neglected.
The electronic structure and magnetic properties of Mott-Hubbard systems are known to depend sensitively on the choice of the Coulomb interaction [33,34]. We therefore check different Hubbard U values, to make our DFT+DMFT based conclusions robust. In particular, we performed DFT+DMFT calculations with a sufficiently smaller Hubbard parameter U = 5 eV (while using the same Hund's coupling, J = 0.89 eV). The DFT+DMFT approach was implemented within the plane-wave pseudopotentials with the generalized gradient approximation in DFT. We employed the fully localized double-counting correction, evaluated from the self-consistently determined local occupancies, to account for the electronic interactions already described by DFT. The spectral functions were computed using the maximum entropy method. Further technical details about the method used can be found in Ref. [25,26].